options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

stan basic

normal distribution

ex1-0.stan

normal distribution

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
mdl=cmdstan_model('./ex1-0.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.16    1.48    2.25    2.21    2.93    4.20
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -25.8051 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       -6.3269   0.000301942   0.000252394      0.9807      0.9807       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##      lp__    -6.33
##      m        2.21
##      s        1.14
##      y1       2.91
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
mcmc
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -7.31  -6.97 1.13 0.83 -9.41 -6.23 1.00      593      596
##      m     2.23   2.22 0.48 0.43  1.48  3.01 1.00      627      642
##      s     1.43   1.34 0.42 0.35  0.93  2.22 1.00     1127      990
##      y1    2.24   2.23 1.57 1.40 -0.30  4.77 1.00     1676     1643
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"    "y1"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"    "y1"
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -7.31  -6.97 1.13 0.83 -9.41 -6.23 1.00      593      596
##      m     2.23   2.22 0.48 0.43  1.48  3.01 1.00      627      642
##      s     1.43   1.34 0.42 0.35  0.93  2.22 1.00     1127      990
##      y1    2.24   2.23 1.57 1.40 -0.30  4.77 1.00     1676     1643

poisson distribution

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    2.00    2.50    2.85    4.00    6.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

in case fit poisson dist. to normal dist.

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -174.991 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -16.5622   2.48557e-05   0.000216293           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -16.56
##      m        2.85
##      s        1.39
##      y1       1.96
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -16.22 -15.91 1.05 0.75 -18.41 -15.22 1.00      882      990
##      m      2.84   2.84 0.34 0.33   2.28   3.40 1.00     1580     1197
##      s      1.52   1.48 0.27 0.24   1.15   2.04 1.00      851      860
##      y1     2.80   2.80 1.57 1.57   0.21   5.29 1.00     2134     2016

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

fit to poisson dist.

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  int y1;
  y1=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -18.63 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       2.69718    0.00015116   4.20026e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     2.70
##      l        2.85
##      y1       0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 3.25   3.52 0.71 0.32 1.88 3.75 1.00     1039     1194
##      l    2.91   2.90 0.38 0.38 2.28 3.56 1.00      735     1050
##      y1   2.84   3.00 1.71 1.48 0.00 6.00 1.00     1833     1969

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

difference test

ex3.stan

mean difference of 2 normal distributions

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -13861 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -35.6979   0.000837279   0.000622165           1           1       48    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -35.70
##      ma       9.78
##      sa       1.12
##      mb      10.58
##      sb       1.96
##      d        0.80
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -36.98 -36.65 1.48 1.30 -39.68 -35.25 1.00      802      970
##      ma     9.78   9.79 0.27 0.27   9.35  10.21 1.00     1611     1325
##      sa     1.22   1.20 0.22 0.21   0.93   1.62 1.00     2402     1684
##      mb    10.58  10.59 0.49 0.48   9.78  11.36 1.00      616      489
##      sb     2.16   2.10 0.40 0.35   1.66   2.90 1.00     1729     1391
##      d      0.80   0.80 0.56 0.53  -0.15   1.72 1.00      698      814

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration     1    2    3    4
##         1  0.66 0.53 0.42 0.49
##         2 -0.25 0.63 1.33 0.93
##         3 -0.36 0.60 1.14 0.83
##         4  0.68 0.66 0.90 0.53
##         5  0.45 0.67 1.32 1.07
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.916

normal regression

ex4-1.stan

single normal regression

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -4.34345e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSt0wv7/model-a38248be373.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       29      -50.2516   0.000701789   0.000697943           1           1       84    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -50.25
##      b0      10.52
##      b1       3.04
##      s        7.48
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -49.89 -49.49 1.45 1.13 -52.79 -48.42 1.01      310      306
##      b0    10.32  10.38 3.86 3.41   3.71  16.62 1.04      105       85
##      b1     3.04   3.04 0.08 0.07   2.92   3.16 1.03      139      132
##      s      8.50   8.26 1.65 1.45   6.27  11.53 1.00      714      459

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

single normal regression, prediction from new explanatory variable x1

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.89582e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       49      -50.2516   0.000264535   0.000143774      0.1514           1      104    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -50.25
##    b0        10.52
##    b1         3.04
##    s          7.48
##    m[1]      15.50
##    m[2]      46.15
##    m[3]      76.79
##    m[4]     107.44
##    m[5]     138.09
##    m[6]     168.73
##    m[7]     199.38
##    m[8]     230.03
##    m[9]     260.67
##    m[10]    291.32
##    y1[1]     13.75
##    y1[2]     57.26
##    y1[3]     72.47
##    y1[4]     94.83
##    y1[5]    132.97
##    y1[6]    172.74
##    y1[7]    198.18
##    y1[8]    207.57
##    y1[9]    257.93
##    y1[10]   297.49
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -49.85 -49.49 1.33 1.05 -52.42 -48.42 1.01      322      298
##    b0      10.59  10.51 4.03 3.69   4.37  17.69 1.02      214      143
##    b1       3.04   3.04 0.08 0.07   2.91   3.16 1.01      266      224
##    s        8.48   8.30 1.49 1.41   6.42  11.21 1.00      637     1083
##    m[1]    15.57  15.49 3.92 3.61   9.50  22.48 1.02      214      143
##    m[2]    46.20  46.20 3.27 3.04  41.10  51.87 1.02      216      134
##    m[3]    76.83  76.87 2.69 2.54  72.64  81.46 1.02      228      185
##    m[4]   107.46 107.48 2.23 2.14 104.00 111.04 1.01      277      340
##    m[5]   138.10 138.09 1.97 1.93 134.89 141.37 1.01      491      881
##    m[6]   168.73 168.71 2.00 1.97 165.42 172.01 1.00     1573     1313
##    m[7]   199.36 199.31 2.30 2.34 195.56 203.15 1.00     2253     1537
##    m[8]   230.00 230.02 2.79 2.76 225.35 234.58 1.00     1479     1243
##    m[9]   260.63 260.60 3.39 3.32 255.17 266.10 1.00      885      740
##    m[10]  291.26 291.27 4.05 3.96 284.87 297.78 1.00      615      604
##    y1[1]   15.74  15.77 9.56 9.21   0.46  31.26 1.00      740     1342
##    y1[2]   46.14  45.89 9.49 9.02  30.69  61.30 1.01      944     1562
##    y1[3]   76.64  76.66 9.29 8.72  61.64  92.04 1.00     1463     1702
##    y1[4]  107.67 107.82 8.80 8.62  93.34 121.44 1.00     1608     1785
##    y1[5]  137.99 138.15 8.85 8.20 123.56 151.94 1.00     1622     1737
##    y1[6]  168.59 168.67 8.92 8.58 153.36 183.19 1.00     1951     1892
##    y1[7]  198.92 198.80 8.85 8.94 184.93 213.30 1.00     1938     1937
##    y1[8]  230.17 230.19 9.34 9.02 214.18 245.01 1.00     2132     1964
##    y1[9]  260.66 260.79 9.15 8.58 245.69 275.14 1.00     1913     1871
##    y1[10] 291.48 291.27 9.44 9.16 276.06 307.20 1.00     1525     1648

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      15.6   15.5  3.92  3.61   9.50  22.5  1.02     214.     143.
##  2 m[2]      46.2   46.2  3.27  3.04  41.1   51.9  1.02     216.     135.
##  3 m[3]      76.8   76.9  2.69  2.54  72.6   81.5  1.02     229.     186.
##  4 m[4]     107.   107.   2.23  2.14 104.   111.   1.01     278.     340.
##  5 m[5]     138.   138.   1.97  1.93 135.   141.   1.01     492.     882.
##  6 m[6]     169.   169.   2.00  1.97 165.   172.   1.00    1574.    1313.
##  7 m[7]     199.   199.   2.30  2.34 196.   203.   1.00    2253.    1538.
##  8 m[8]     230.   230.   2.79  2.76 225.   235.   1.00    1480.    1243.
##  9 m[9]     261.   261.   3.39  3.32 255.   266.   1.00     886.     741.
## 10 m[10]    291.   291.   4.05  3.96 285.   298.   1.00     615.     604.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     15.7   15.8  9.56  9.21   0.459  31.3 1.00      740.    1342.
##  2 y1[2]     46.1   45.9  9.49  9.02  30.7    61.3 1.01      945.    1563.
##  3 y1[3]     76.6   76.7  9.29  8.72  61.6    92.0 1.00     1463.    1702.
##  4 y1[4]    108.   108.   8.80  8.62  93.3   121.  1.00     1609.    1785.
##  5 y1[5]    138.   138.   8.85  8.20 124.    152.  1.00     1622.    1737.
##  6 y1[6]    169.   169.   8.92  8.58 153.    183.  1.00     1952.    1893.
##  7 y1[7]    199.   199.   8.85  8.94 185.    213.  0.999    1939.    1938.
##  8 y1[8]    230.   230.   9.34  9.02 214.    245.  1.00     2133.    1964.
##  9 y1[9]    261.   261.   9.15  8.58 246.    275.  1.00     1914.    1871.
## 10 y1[10]   291.   291.   9.44  9.16 276.    307.  1.00     1525.    1649.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

single normal regression, prediction from original explanatory variable x

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -49.84 -49.51 1.29 1.08 -52.42 -48.42 1.00      495      534
##    b0      10.54  10.66 3.68 3.41   4.63  16.37 1.01      222      223
##    b1       3.04   3.03 0.07 0.07   2.92   3.15 1.01      297      428
##    s        8.47   8.29 1.51 1.37   6.39  11.30 1.00     1030      848
##    m[1]   140.85 140.85 1.99 1.90 137.49 143.98 1.00      685      921
##    m[2]   255.85 255.84 3.35 3.28 250.01 261.11 1.00     1072     1166
##    m[3]    16.72  16.81 3.56 3.30  10.97  22.33 1.01      222      225
##    m[4]    56.66  56.75 2.82 2.64  52.11  61.17 1.01      233      238
##    m[5]   218.20 218.24 2.69 2.67 213.64 222.41 1.00     2060     1440
##    m[6]    15.51  15.61 3.58 3.31   9.75  21.15 1.01      222      225
##    m[7]   213.08 213.15 2.61 2.55 208.65 217.15 1.00     2168     1519
##    m[8]   148.88 148.86 1.99 1.90 145.52 152.00 1.00      896     1146
##    m[9]    83.91  84.02 2.41 2.29  79.98  87.80 1.01      264      313
##    m[10]   97.86  97.92 2.24 2.13  94.18 101.39 1.01      296      433
##    m[11]  236.23 236.24 2.99 2.97 230.98 240.92 1.00     1460     1311
##    m[12]   72.22  72.34 2.58 2.43  67.98  76.37 1.01      247      256
##    m[13]   93.37  93.43 2.29 2.17  89.62  97.01 1.01      284      385
##    m[14]   88.86  88.93 2.34 2.23  85.04  92.59 1.01      273      332
##    m[15]  103.85 103.91 2.18 2.10 100.26 107.25 1.01      316      443
##    m[16]  210.47 210.53 2.57 2.50 206.12 214.48 1.00     2222     1520
##    m[17]  291.17 291.18 4.05 3.95 284.29 297.59 1.00      754      949
##    m[18]  203.46 203.50 2.47 2.40 199.30 207.34 1.00     2304     1496
##    m[19]  130.83 130.81 2.00 1.89 127.45 133.94 1.00      518      756
##    m[20]  229.93 229.96 2.88 2.86 224.98 234.45 1.00     1646     1370
##    y1[1]  140.87 140.80 8.59 8.21 127.16 155.12 1.00     2031     1650
##    y1[2]  255.86 255.66 9.48 9.19 240.35 271.57 1.00     1683     1732
##    y1[3]   16.66  16.90 9.24 9.05   1.46  31.96 1.00     1080     1446
##    y1[4]   56.75  56.84 9.22 8.88  42.54  72.00 1.00     1223     1788
##    y1[5]  218.14 218.27 9.09 8.69 203.36 232.98 1.00     2180     1884
##    y1[6]   15.22  15.52 9.31 8.86  -0.80  30.36 1.00      876     1357
##    y1[7]  213.51 213.57 9.04 8.85 198.64 228.31 1.00     2142     1758
##    y1[8]  148.90 148.88 8.92 8.48 133.98 163.38 1.00     1877     1869
##    y1[9]   84.08  83.81 8.96 8.62  69.58  98.84 1.00     1555     1711
##    y1[10]  97.75  97.65 8.90 8.85  83.16 112.47 1.00     1476     1645
##    y1[11] 236.10 236.07 9.09 8.86 221.04 250.92 1.00     1678     1493
##    y1[12]  71.99  72.23 9.05 8.63  57.11  86.79 1.00     1311     1522
##    y1[13]  93.28  93.32 9.05 8.48  78.02 107.96 1.00     1558     1622
##    y1[14]  88.88  88.96 8.87 8.33  74.14 103.85 1.00     1691     1813
##    y1[15] 103.96 103.72 8.85 8.45  89.81 118.65 1.00     1886     1935
##    y1[16] 210.08 210.37 8.79 8.49 195.21 224.62 1.00     1831     1749
##    y1[17] 290.90 291.08 9.57 9.11 275.56 306.53 1.00     1760     1735
##    y1[18] 203.51 203.47 9.01 8.67 188.48 217.93 1.00     1626     1892
##    y1[19] 131.12 130.96 8.64 8.12 116.99 145.29 1.00     1539     1912
##    y1[20] 229.91 230.18 8.99 8.41 215.21 244.50 1.00     1986     1892

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation

ex4-41.stan

use covariance matrix and multinormal distribution

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

use covariance matrix and wishart distribution

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

use correlation matrix and wishart distribution

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan use covariance matrix and multinormal distribution

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##       [,1]  [,2]
## [1,] 0.857 0.664
## [2,] 0.664 1.158
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.666
## [2,] 0.666 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -120.143 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -13.0271   3.32791e-05   5.95564e-05           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -13.03
##   m[1]       -0.18
##   m[2]       -0.18
##   s[1]        0.90
##   s[2]        1.05
##   r           0.67
##   cv[1,1]     0.81
##   cv[2,1]     0.63
##   cv[1,2]     0.63
##   cv[2,2]     1.10
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -17.13 -16.73 1.81 1.64 -20.35 -14.95 1.01      699     1006
##   m[1]     -0.18  -0.19 0.24 0.23  -0.57   0.21 1.00     1137     1125
##   m[2]     -0.19  -0.20 0.28 0.26  -0.62   0.27 1.00     1101      968
##   s[1]      1.01   0.98 0.18 0.17   0.75   1.34 1.00     1372     1342
##   s[2]      1.17   1.13 0.21 0.19   0.88   1.56 1.00     1311     1465
##   r         0.62   0.64 0.15 0.14   0.34   0.82 1.00      614      831
##   cv[1,1]   1.05   0.96 0.41 0.32   0.56   1.79 1.00     1372     1342
##   cv[2,1]   0.76   0.69 0.37 0.30   0.30   1.44 1.00      711     1031
##   cv[1,2]   0.76   0.69 0.37 0.30   0.30   1.44 1.00      711     1031
##   cv[2,2]   1.41   1.29 0.55 0.42   0.78   2.43 1.00     1311     1465

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -52.4202 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -13.3503   0.000292043   3.38327e-05           1           1       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -13.35
##   s[1]        0.93
##   s[2]        1.08
##   r           0.67
##   cv[1,1]     0.86
##   cv[2,1]     0.66
##   cv[1,2]     0.66
##   cv[2,2]     1.16
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -16.17 -15.81 1.32 1.01 -18.79 -14.79 1.00      874      978
##   s[1]      1.01   0.98 0.19 0.17   0.75   1.33 1.00     1088     1036
##   s[2]      1.17   1.14 0.22 0.19   0.88   1.55 1.00     1114     1260
##   r         0.62   0.64 0.14 0.14   0.37   0.82 1.00      500      693
##   cv[1,1]   1.05   0.96 0.41 0.32   0.56   1.77 1.00     1088     1036
##   cv[2,1]   0.76   0.70 0.38 0.30   0.33   1.49 1.00      581      793
##   cv[1,2]   0.76   0.70 0.38 0.30   0.33   1.49 1.00      581      793
##   cv[2,2]   1.41   1.29 0.57 0.42   0.78   2.41 1.00     1114     1260

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -36.1841 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -13.3503   0.000599028   0.000839941           1           1       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -13.35
##   s[1]        0.93
##   s[2]        1.08
##   r[1,1]      1.00
##   r[2,1]      0.67
##   r[1,2]      0.67
##   r[2,2]      1.00
##   cv[1,1]     0.86
##   cv[2,1]     0.66
##   cv[1,2]     0.66
##   cv[2,2]     1.16
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -15.57 -15.23 1.36 1.10 -18.32 -14.10 1.00      775     1110
##   s[1]      1.01   0.98 0.19 0.18   0.76   1.33 1.00     1220      906
##   s[2]      1.17   1.15 0.22 0.20   0.88   1.57 1.00      946     1030
##   r[1,1]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   r[2,1]    0.61   0.64 0.15 0.14   0.32   0.83 1.00      741      891
##   r[1,2]    0.61   0.64 0.15 0.14   0.32   0.83 1.00      741      891
##   r[2,2]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   cv[1,1]   1.05   0.97 0.43 0.34   0.57   1.78 1.00     1220      906
##   cv[2,1]   0.76   0.68 0.41 0.32   0.29   1.52 1.00      719      663
##   cv[1,2]   0.76   0.68 0.41 0.32   0.29   1.52 1.00      719      663
##   cv[2,2]   1.43   1.32 0.58 0.45   0.77   2.47 1.00      946     1030

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -624.153 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -13.0271   0.000339274   0.000683859           1           1       16    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__       -13.03
##  m[1]        -0.18
##  m[2]        -0.18
##  cov[1,1]     0.81
##  cov[2,1]     0.63
##  cov[1,2]     0.63
##  cov[2,2]     1.10
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -15.00 -14.65 1.74 1.58 -18.35 -12.85 1.00      761      933
##  m[1]      -0.18  -0.18 0.25 0.24  -0.60   0.23 1.00     1338     1237
##  m[2]      -0.18  -0.18 0.29 0.28  -0.64   0.28 1.00     1162     1339
##  cov[1,1]   1.25   1.12 0.56 0.40   0.66   2.25 1.00     1444     1098
##  cov[2,1]   0.97   0.86 0.53 0.37   0.35   1.93 1.00     1085      817
##  cov[1,2]   0.97   0.86 0.53 0.37   0.35   1.93 1.00     1085      817
##  cov[2,2]   1.67   1.50 0.70 0.49   0.90   3.00 1.00     1223      850

distributions

normal dist.

distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2

ex1-0.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
n=20
y=rnorm(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.53    1.45    1.98    1.98    2.47    3.21
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -60.248 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13       -3.5906   3.72627e-05   4.97546e-05           1           1       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -3.59
##      m        1.98
##      s        0.73
##      y1       2.96
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -4.95  -4.60 1.11 0.70 -7.22 -3.95 1.00      762      920
##      m     1.98   1.99 0.18 0.17  1.70  2.26 1.00     1628     1082
##      s     0.80   0.78 0.14 0.13  0.61  1.07 1.00     1682     1166
##      y1    1.96   1.97 0.86 0.84  0.58  3.34 1.00     1999     1932

Bernoulli dist.

The event with probablity p occur y=1, not occur y=0 
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)

ex1-1.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  y~bernoulli(p);
}
generated quantities{
  real s;
  s=(p*(1-p))^.5;
}
n=10
y=sample(0:1,n,replace=T)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00    0.30    0.75    1.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -7.53655 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -6.10864    0.00573752   0.000129064           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -6.11
##      p        0.30
##      s        0.46
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -8.13  -7.86 0.68 0.31 -9.53 -7.64 1.00      916      957
##      p     0.33   0.32 0.13 0.13  0.14  0.55 1.00      729      740
##      s     0.45   0.47 0.05 0.04  0.34  0.50 1.00      719      740

geometric dist.

The event with probablity p occur after y trials(failure) 
y~Ge(p), y[0,Inf), P(Y=y)=p(1−p)^y
E[y]=1/p, V[y]=(1-p)/p^2

ex1-2.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  for (n in 1:N) {
    target+=log(p)+y[n]*log(1-p);
  }
}
generated quantities{
  real m;
  m=1/p;
  real s;
  s=((1-p)/p^2)^.5;
}
n=20
y=rgeom(n,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.75    1.00    2.20    4.00    8.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -45.6597 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -39.7495       0.00348   0.000336829           1           1        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -39.75
##      p        0.31
##      m        3.20
##      s        2.65
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -41.80 -41.51 0.72 0.31 -43.22 -41.28 1.00     1117     1403
##      p      0.32   0.32 0.06 0.06   0.23   0.42 1.00      759     1010
##      m      3.22   3.13 0.61 0.55   2.36   4.39 1.00      759     1010
##      s      2.67   2.58 0.62 0.56   1.80   3.85 1.00      759     1010

negative binomial dist.

The event with probablity p occur a times after y-1 trials 
y~NB(a,p), y[a,Inf), P(Y=y)=p^a*(1−p)^(y-a-1)
E[y]=a*(1-p)/p, V[y]=a*(1-p)/p^2

ex1-3.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
  real<lower=1> a;
}
model{
  y~neg_binomial(a,p);
}
generated quantities{
  real m;
  m=a*(1-p)/p;
  real s;
  s=(a*(1-p)/p^2)^.5;
  int y1;
  y1=neg_binomial_rng(a,p);
}
n=50
a=3
y=rnbinom(n,a,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    4.00    6.50    7.38    9.00   23.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -150.761 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13       -142.64    0.00340739    0.00395503      0.9713      0.9713       16    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -142.64
##      p        0.52
##      a        3.81
##      m        3.57
##      s        2.63
##      y1       7.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -144.10 -143.77 1.13 0.84 -146.39 -142.99 1.01      389      775
##      p       0.63    0.61 0.17 0.20    0.37    0.93 1.01      248      453
##      a       4.64    4.53 1.25 1.45    2.79    6.82 1.01      255      447
##      m       2.75    2.82 1.37 1.54    0.50    4.89 1.01      264      455
##      s       2.15    2.15 0.89 0.96    0.74    3.60 1.01      254      455
##      y1      7.45    7.00 4.47 4.45    1.00   15.00 1.00     2076     1647

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

binomial dist.

The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)

ex2-3.stan

data{
  int N;
  int n;
  array[N] int y;
}
parameters{
  real<lower=0,upper=0> p;
}
model{
  y~binomial(n,p);
}
generated quantities{
  real m;
  m=n*p;
  real s;
  s=(n*p*(1-p))^.5;
  int y1;
  y1=binomial_rng(n,p);
}
n=20
n0=10
y=rbinom(n,n0,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    2.00    3.00    2.95    4.00    5.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -136.737 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -121.314    0.00270906   0.000861894           1           1        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -121.31
##      p        0.30
##      m        2.95
##      s        1.44
##      y1       2.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -123.34 -123.11 0.63 0.30 -124.53 -122.88 1.00     1130     1339
##      p       0.30    0.30 0.03 0.03    0.25    0.35 1.00      879     1122
##      m       2.98    2.97 0.31 0.32    2.49    3.49 1.00      879     1122
##      s       1.44    1.45 0.04 0.04    1.37    1.51 1.00      879     1122
##      y1      2.95    3.00 1.44 1.48    1.00    6.00 1.00     1763     1701

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

ex2-31.stan

trial n is varied from each sample

data{
  int N;
  array[N] int n;
  array[N] int y;
}
parameters{
  real<lower=0> p;
}
model{
  y~binomial(n,p);
}
n=20
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-31.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -49.9164 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        3      -49.3093    0.00401524   0.000319669      0.9516      0.9516        5    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -49.31
##      p        0.32
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -51.35 -51.09 0.69 0.35 -52.78 -50.84 1.00     1144     1428
##      p      0.32   0.32 0.05 0.05   0.24   0.41 1.00      809     1124

Poisson dist.

The event occur l times in unit range, the event occur y times. 
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=l

ex1-4.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0> l;
}
model{
  y~poisson(l);
}
generated quantities{
  real s;
  s=l^.5;
  int y1;
  y1=poisson_rng(l);
}
n=20
y=rpois(n,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     2.0     3.0     3.2     4.0     7.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-4.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -22.8192 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       10.4417   0.000665407   0.000481145           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    10.44
##      l        3.20
##      s        1.79
##      y1       2.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad   q5   q95 rhat ess_bulk ess_tail
##      lp__ 11.15  11.40 0.63 0.29 9.84 11.61 1.01     1012     1316
##      l     3.27   3.26 0.39 0.39 2.65  3.93 1.00      803     1224
##      s     1.81   1.81 0.11 0.11 1.63  1.98 1.00      803     1224
##      y1    3.28   3.00 1.86 1.48 1.00  7.00 1.00     1888     1944

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

exponential dist.

The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2

ex1-5.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> l;
}
model{
  y~exponential(l);
}
generated quantities{
  real m;
  m=1/l;
  real s;
  s=1/l;
  real y1;
  y1=exponential_rng(l);
}
n=20
y=rexp(n,2)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.007   0.095   0.378   0.448   0.537   1.920
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -34.0293 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -3.94689   0.000679157   1.92585e-05           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -3.95
##      l        2.23
##      m        0.45
##      s        0.45
##      y1       0.44
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -3.62  -3.34 0.71 0.30 -5.07 -3.12 1.00      846     1186
##      l     2.35   2.30 0.51 0.49  1.57  3.27 1.00      741     1017
##      m     0.45   0.44 0.10 0.09  0.31  0.64 1.00      741     1017
##      s     0.45   0.44 0.10 0.09  0.31  0.64 1.00      741     1017
##      y1    0.45   0.30 0.47 0.32  0.02  1.36 1.00     1853     1860

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

gamma dist.

The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2

ex1-6.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> a;
  real<lower=0> l;
}
model{
  y~gamma(a,l);
}
n=20
y=rgamma(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.41    0.79    1.60    2.03    2.55    6.28
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -141.353 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12      -32.3719   0.000120154   5.75695e-06           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -32.37
##      a        1.83
##      l        0.90
##      m        2.03
##      s        1.50
##      y1       1.10
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -32.73 -32.45 0.98 0.78 -34.69 -31.75 1.01      669     1074
##      a      2.12   2.07 0.59 0.59   1.23   3.10 1.01      617      728
##      l      1.07   1.04 0.33 0.32   0.57   1.64 1.01      618      511
##      m      2.03   2.00 0.33 0.32   1.56   2.61 1.00     2163     1619
##      s      1.43   1.38 0.32 0.27   1.03   2.05 1.01      794      844
##      y1     2.03   1.71 1.46 1.22   0.33   4.77 1.00     1849     1735

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

log normal dist.

logalithm of variable follows normal distribution
log y~N(m,s), y>0
E[y]=exp(m+s^2)
V[y]=exp(2*m+s^2)*(exp(s^2)-1)

ex1-7.stan

data{
  int N;
  vector[N]<lower=0> y;
}
parameters{
  real m0;
  real<lower=0> s0;
}
model{
  y~lognormal(m0,s0);
}
n=30
y=rlnorm(n,0.5,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.2     1.2     1.6     3.2     3.1    35.6
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -55.1291 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -40.1814    0.00103138    0.00017249           1           1       11    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -40.18
##      m0       0.56
##      s0       0.92
##      m        4.11
##      s        3.12
##      y1       2.01
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -41.29 -41.01 1.02 0.77 -43.26 -40.31 1.00      987     1058
##      m0     0.57   0.57 0.19 0.18   0.26   0.86 1.00     1282     1251
##      s0     0.99   0.97 0.13 0.13   0.79   1.22 1.00     1894     1316
##      m      5.05   4.54 1.99 1.38   2.97   8.70 1.00     1600     1334
##      s      4.05   3.55 1.93 1.32   2.07   7.53 1.00     1668     1383
##      y1     2.97   1.80 3.82 1.58   0.34   9.46 1.00     1945     2054

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

beta dist.

use as prior of binomial dist parameter p
event with probabilty p0 occur a-1 times, do not occur b-1 times in a+b-2 trials
p~Be(a,b), p[0,1], p=p0^(a-1)*(1-p0)^(b-1)
E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)

ex1-8.stan

data{
  int N;
  vector[N] p;
}
parameters{
  real<lower=0> a;
  real<lower=0> b;
}
model{
  p~beta(a,b);
}
n=20
p=rbeta(n,1,2)
summary(p)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.006   0.225   0.293   0.301   0.339   0.607
par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')

data=list(N=n,p=p)

mdl=cmdstan_model('./ex1-8.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -17.3024 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7       9.06616    0.00223233   4.62945e-05      0.9944      0.9944       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     9.07
##      a        1.97
##      b        4.71
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad   q5   q95 rhat ess_bulk ess_tail
##      lp__ 10.48  10.76 0.95 0.69 8.46 11.40 1.00      565      803
##      a     2.23   2.16 0.61 0.60 1.32  3.31 1.00      652      748
##      b     5.39   5.23 1.55 1.57 3.11  8.20 1.01      603      657

Dirichlet dist.

use as prior of categorical dist parameter p
p~dir(p0), p[0,1]

ex1-9.stan

data {
  int N;
  int K;
  matrix[N, K] p;
}
parameters {
  simplex[K] p0;
}
model {
  for(i in 1:N){
     p[i]~dirichlet(p0);
  }
}
library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)
##        V1              V2              V3       
##  Min.   :0.003   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.264   1st Qu.:0.096   1st Qu.:0.001  
##  Median :0.632   Median :0.314   Median :0.008  
##  Mean   :0.530   Mean   :0.331   Mean   :0.139  
##  3rd Qu.:0.762   3rd Qu.:0.483   3rd Qu.:0.117  
##  Max.   :0.962   Max.   :0.983   Max.   :0.923
boxplot(p)

data=list(N=n,K=ncol(p),p=p)

mdl=cmdstan_model('./ex1-9.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = 38.3768 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       56.0457   0.000115464   1.19815e-05           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__     56.05
##     p0[1]     0.48
##     p0[2]     0.34
##     p0[3]     0.18
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__  51.50  51.84 1.07 0.70 49.37 52.47 1.00     1013     1180
##     p0[1]  0.48   0.48 0.06 0.06  0.37  0.57 1.00     1866     1396
##     p0[2]  0.34   0.34 0.06 0.06  0.25  0.43 1.00     1444     1147
##     p0[3]  0.19   0.18 0.04 0.04  0.13  0.25 1.00     1847     1535